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We develop a simple but rigorous model of protein-protein association kinetics based on diffusional 
association on free energy landscapes obtained by sampling configurations within and surrounding 
the native complex binding funnels. Guided by results obtained on exactly solvable model problems, 
we transform the problem of diffusion in a potential into free diffusion in the presence of an absorbing 
zone spanning the entrance to the binding funnel. The free diffusion problem is solved using a 
recently derived analytic expression for the rate of association of asymmetrically oriented molecules. 
Despite the required high steric specificity and the absence of long-range attractive interactions, the 
computed rates are typically on the order of 10*-10® s~^, several orders of magnitude higher 

than rates obtained using a purely probabilistic model in which the association rate for free diffusion 
of uniformly reactive molecules is multiplied by the probability of a correct alignment of the two 
partners in a random collision. As the association rates of many protein-protein complexes are also in 
the lO'^-lO'' s-\ our results suggest that free energy barriers arising from desolvation and/or 

side-chain freezing during complex formation or increased ruggedness within the binding funnel, 
which are completely neglected in our simple diffusional model, do not contribute significantly to the 
dynamics of protein-protein association. The transparent physical interpretation of our approach 
that computes association rates directly from the size and geometry of protein-protein binding 
funnels makes it a useful complement to Brownian dynamics simulations. 
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I. INTRODUCTION 

The calculation of rates of protein-protein association 
is of great interest to biology. These rates span a wide 
range of values, from approximately 10^ to 10^° s~^. 
If the two proteins are modeled as uniformly reactive 
spheres, the diffusion-limited rate constant is simply 
given by the classical Smoluchowski expression (Smolu- 
chowski 1917), fcon = 47rZ3i? (where D is the relative 
translational diffusion constant and R is the sum of the 
radii), which yields rates of 10^-10^° s~^ for asso- 
ciations relevant to proteins. Usually, however, proteins 
exhibit a highly anisotropic distribution of reactivity over 
their surface. This can be modelled by localized reactive 
sites on the surface of the proteins that have to be suf- 
ficiently precisely aligned for the complex formation to 
occur. 

Purely probabilistic models have tried to account for 
such steric constraints by multiplying the Smoluchowski 
rate for uniform spheres by the probability that, in a ran- 
dom encounter, the two molecules are properly aligned 
("geometric rate") (Janin (1997) gives an example of this 
method). This yields rate constants which are typically 
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several orders of magnitude lower than the Smoluchowski 
diffusion-limited rate and are usually much smaller than 
the values experimentally observed for biological com- 
plexes. It has been found that this discrepancy can 
be moderated by taking into account the effect of ro- 
tational diffusion (Shoup et al. 1981, Northrup and Er- 
ickson 1992); additional rate enhancements are brought 
about by the presence of attractive interparticle forces 
("electrostatic steering"; see Schreiber and Fersht 1996, 
GabdouUine and Wade 1997, Vijayakumar et al. 1998), 
and the formation of a weakly specific, loosely bound en- 
counter complex that subsequently evolves into the final 
bound state (Selzer and Schreiber 1999, Camacho et al. 
2000). 

To replace the estimation of protein-protein associ- 
ation rates via the geometric rate by a more accurate 
method, most authors have pursued a computational ap- 
proach by carrying out explicit numerical simulations of 
the diffusional association of macromolecules, commonly 
referred to as Brownian dynamics (BD) simulations (for 
an excellent review, see GabdouUine and Wade 2002). 
Here, the protein molecules are modeled in varying de- 
tail, from a simple spherical approximation up to full 
atomic detail. In the simulation, the molecules are ini- 
tially placed in random orientations at a fixed initial 
separation b. Diffusional trajectories, with or without 
the presence of an interparticle force (such as electro- 
static interactions), are then generated by means of the 
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Ermak-McCammon algorithm (Ermak and McCammon 
1978). A trajectory is ended either when the molecules 
have come together in proper orientation to successfully 
form a complex, or when their separation has exceeded 
a certain truncation value c > b such that the probably 
for an encounter has become vanishingly small. The frac- 
tion of 'successful' trajectories is then used to compute 
the association rate fcon- 

Northrup and Erickson (1992) have used such a BD 
simulation to compute the association rate of spherical 
molecules with a reactive patch, consisting of four con- 
tact points in a 17A X 17A square arrangement on a 
plane tangential to the surface of the molecules. Reac- 
tion is then assumed to occur if three of the four con- 
tact points are correctly matched and within a specified 
maximum distance. In the absence of any interparticle 
forces, the authors find an association rate of fcon — 10^ 
M"-'^ s~^, about two orders of magnitude higher than 
the geometric rate. GabdouUine and Wade (2001) com- 
pute association rates for five protein-protein complexes 
using full-atom structures in the presence of long-range 
electrostatic forces. The reaction condition is defined by 
formation of subsets of the polar contacts observed in the 
native complex structure. 

In this paper, we present a different route towards es- 
timation of rates of bimolecular association. Instead of 
employing a computer simulation based approach such 
as the method of BD simulations outlined above, we 
use a recently derived analytical expression (Schlosshauer 
and Baker 2002) for the association rate of two spherical 
molecules with anisotropic reactivity in the absence of 
any interaction forces. The reaction condition is formu- 
lated by specifying the ranges of mutual orientations of 
the two molecules for which complex formation will oc- 
cur. We thus do not require an exact mutual alignment 
of the binding partners, but instead assume that favor- 
able short-range interactions "guide" the molecules into 
their final bound configurations once the molecules are 
oriented within specified angular tolerances (see Fig.^). 
These tolerances can therefore be viewed as an implicit 
modelling of attractive short-range forces. We derive es- 
timates for the tolerances from free energy landscapes ob- 
tained by sampling configurations within and surround- 
ing the native binding funnel. These values arc then used 
in our analytical expression to compute the correspond- 
ing association rates. By determining the size and geom- 
etry of the aperture in phase space which must be entered 
for binding to occur, and rigorously solving the problem 
of diffusion through this aperture, our approach provides 
a physically transparent complement to BD simulations 
for computing binding rates from structures of protein- 
protein complexes. 



II. RESULTS 

To compute protein-protein association rates from 
the three-dimensional structures of protein-protein com- 




FIG. 1: Simple model of binding dynamics. Attractive short- 
range forces produce a funnel in the free energy landscape 
leading into the native complex. Once the molecules descend 
several kT into the funnel, they are effectively captured and 
binding occurs rapidly. In our simple model, the rate of as- 
sociation is approximated by the rate of free diffusion into 
a reactive zone in phase space, as indicated schematically in 
the XY plane of the drawing. To compute the rate of associa- 
tion, we need to first determine the dimensions of the reactive 
zone, and second, compute the rate of free diffusion into this 
zone. A more general model would include long-range (elec- 
trostatic) interactions which would bias the diffusion process 
towards the funnel entrance. 



plexes according to the simple diffusive model described 
above and in Fig. 1, three ingredients are required. The 
first is a general theory for computing the diffusion- 
limited association rate as a function of the orientational 
constraints associated with properly aligigning the two 
binding sites (the size and shape of the reactive zone in 
Fig. 1). The second is a method for transforming a diffu- 
sion in a potential problem into a free diffusion problem — 
in the context of Fig. 1, an estimate of how deeply the 
reactive zone lies within the binding funnel (i.e. how far 
molecules must descend into the binding funnel before 
they are effectively captured). The third is a method for 
mapping the binding funnel for two proteins given the 
crystal structure of the protein-protein complex. We ad- 
dress these issues in Sees. I, II, and III, and in Sec. IV 
use the results to compute approximate diffusion limited 
association rates from the structures of protein-protein 
complexes. 



3 




FIG. 2: The axes and angles relevant to the reaction con- 
dition, Eqs. (0. The angles 6a and Ob measure how close 
the center of each reactive patch (coninciding with the re- 
spective body-fixed z axis) is to the center-to-center vector 
(dashed line). The angles Sij) and 5x denote relative torsion 
angles of the two body-fixed coordinate systems {xA,yA, za) 
and (xB,yB, zb)- For the sake of easier visualization of these 
two angles, the origin of the xa and yA axes (belonging to 
the coordinate system of sphere A) has been shifted such as 
to coincide with the origin of the coordinate system of sphere 
B. Our reaction condition, Eqs. iQ, requires near-optimal 
alignment, i.e. all angles 6a, 6b, S(f>, and 5x must be below 
given limits. 



A. Theory for the difFusion-limited association rate 
with general orientational constraints 

Here, we shall restrict ourselves to a brief outline; the 
full derivation of our expression for the association rate 
constant in the presence of general orientational con- 
straints can be found in Schlosshauer and Baker (2002). 

We consider translational and rotational diffusional 
motion of two spherical molecules A and B with radii 
Ra and Rb- To derive an expression for the association 
rate constant, we solve the steady-state translational- 
rotational diffusion equation describing the diffusional 
motion of the two spheres, subject to a reaction condition 
that ensures that binding can only occur if the mutual 
orientation of the two spheres is sufficiently close to the 
orientation in the bound configuration that defines the 
optimal alignment. 

The reaction conditionis implemented as follows (see 
Fig. [5J): The centers of "reactive patches" on the two 
spheres are defined by the intersection of the center-to- 
center vector with the surfaces of the spheres in the native 
bound configuration. Each sphere carries its own body- 
fixed coordinate system {xg, Us, Zg}, s — A, B, where the 
Zs axis points at the center of the reactive patch. The 
angles 0a and 9b then quantify the distance of the cen- 



ter of each reactive patch to the center-to-center vec- 
tor, whereas the angles 6(j) and Sx denote relative tor- 
sional angles between the body-fixed coordinate systems 
{xA,yA,ZA) and {xb,Vb,zb). 

At a first glance, one might assume that a fifth param- 
eter is required to fully describe the mutual orientation 
of the two spheres — namely, an azimuthal angle <f)A in 
addition to the polar angle Oa to fix the location of the 
reactive patch on the surface of sphere A. For the for- 
mulation of the reaction condition, however, four angles 
suffice, because the position of the center of the reactive 
patch is automatically specified through its coincidence 
with the ZA axis. This leaves only one free parameter, 
namely the "width" of the patch, which is described by 
the angle 9a- 

The optimal alignment is then defined hy 6a ^ Ob — 
Scf) = 6x ^ (additionally, the length r of the center- 
to-center vector must be equal to the sum of the radii 
of the spheres). Our reaction condition requires that all 
these angles are suffiently close to zero for the reaction 
to occur, i.e. that the following conditions are fulfilled: 



Ra + Rb = R 



' < S6n 



(1) 



Sx < Sxo 



Using the constant-fiux approximation introduced by 
Shoup et al. (1991), we obtain for the association rate 
constant (Schlosshauer and Baker 2002): 
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diffusion constant, oq — (47r)'^(5(/)o<5xo(l — cofi9\){\ — 
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quantifies the extent of diffusion control in the reaction. 
Furthermore, 
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where dl^j^{9) denotes the Wigner rotation function. 
( L mA ms ) is the Wigner 3-j symbol, and C = 

R[{D'^f/D)lA{lA + 1) + {D'S'/D)Ib{Ib + l)Y'\ where 
Z?™* and ZJg"* are the rotational diffusion constants. 
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flo = ao/(47r x Stt^ x Stt^) represents the fraction of angu- 
lar orientational space over which the reaction can occur, 
and the geometric rate is thus given by fcon = 4:nDR x ao- 



B. Transformation of the diffusion in a potential 
problem into a free diffusion problem 



A crucial point in the application of Eq. (jSJ is the esti- 
mation of the angular constraints 6^, Og, d<po, and Sxo- 
We would like to estimate the ranges in mutual orien- 
tation of the two proteins for which short-range attrac- 
tive forces between the atoms are sufficiently dominant 
to guide the two molecules into the final bound configura- 
tion, and then translate the problem of diffusional asso- 
ciation in the attractive potential into free diffusion with 
an absorbing region in configurational space. To moti- 
vate this mapping, we shall first study two simple toy 
models for translational and rotational diffusion, respec- 
tively. We then use these ideas to explicitly obtain the 
angular constraints for real protein-protein complexes. 

a. Toy model for translational diffusion. The reac- 
tion rate for diffusion-controlled bimolecular association 
of uniformly reactive spheres in the presence of a poten- 
tial U{r) can be calculated exactly and is given by the 
expression 



dre^^^-^/r' 



(4) 



where /3 = and Ri is interpreted as the center-to- 
center distance between the associating partners at which 
the reaction is assumed to occur. Eq. is the classical 
result derived by Debye (1942). In the absence of any 
potential {U = 0), this simplifies to the Smulochowski 
rate constant for free diffusion with an absorbing region 
of width Rq , given by = AttDRq (in the following we 
use the label "0" to refer to the free diffusion problem, 
and the label "1" to refer to the diffusion in a potential 
problem) . 

It is clear that for Rq — kin > fcon since the 
presence of the (attractive) potential will increase the 
association rate. To find the free diffusion analogue of 
diffusional association in an attractive potential, we in- 
crease i?o until /con = fcon • In other words, for a given 
potential, we can determine the size of the absorbing re- 
gion (the "capture radius") required in the case of free 
diffusion to obtain an association rate equivalent to that 
of diffusion in the potential. A similar redefinition of the 
effective absorbing radius to account for the presence of 
the potential was first introduced by Debye (1942). 

A concrete illustration of this procedure is described in 
Fig. O for a Lennard-Jones potential. For a broad range 
of parameter values, we find that it is sufficient to drop 
down by an energy amount of only 0{AE) = kT to enter 
the capture zone, corresponding to less than 5% of the 
total depth of the potential. The capture radius i?o is 
found to be relatively insensitive to the depth e of the 




FIG. 3: Mapping of the problem of diffusion in a po- 
tential onto that of free diffusion with an absorbing region 
for rotational diffusion on a spherical surface. We use an 
attractive Gaussian potential U{9) — — eexp [-{a6f] with 
e = 10 kcal/mol and a = tt, where the latter corresponds to a 
(half) width of the potential of y ln(2)/cr — 15°, a reasonable 
assumption for a short-range potential. Equating the result- 
ing association rate, Eq. 0, with the rate for free diffusion in 
presence of an absorbing region at 9 = 6o, Eq. we obtain 
9o ~ 33° for the width of the absorbing region, corresponding 
to an energy drop of AE = U{tv) — U{9q) ^ 0.4 kcal/mol. 



potential well, whereas its dependence on the width a 
is much stronger — as it must be, since Rq is an indirect 
measure of the range of the potential. 

The model calculations show that the effect of an at- 
tractive potential U (r) on the association rate can be 
effectively represented by an increase in the radius of 
the interacting spheres, but that in the relevant case of 
protein-protein interactions the relative increase is very 
small (about 7% in our example). Since the association 
rate, Eq. |j2Jl, is largely insensitive to small changes in the 
value of R, we conclude that the approximation of using 
a fixed value r = Ra + Rb = R for the center-to-center 
distance of the two proteins required for the reaction to 
occur (see our reaction condition, Eqs. CQl), rather than 
employing a range of allowed values (such as demanding 
that r < Ra + Rb + SR in Eqs. (QJ), is justified. 

b. Toy model for rotational diffusion. As another 
important illustration of our mapping procedure we shall 
consider two-dimensional rotational diffusion on a spher- 
ical surface in an attractive Gaussian potential U{9) = 
— e exp [-{a0f], with e > 0. A gain, we would like to 
translate this problem into that of free rotational diffu- 
sion with an absorbing region at 9 ~ 9o. 

Solving the rotational diffusion equation in presence of 
a potential U{9) yields in the diffusion-controlled limit 



(5) 



where we let 6*1 ^ for [/ ^0 (diffusion in potential). 
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~ Gaussian potential U(e)^ - e exp[-(aO) ] 



, 9q (free diffusion) 
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FIG. 4: Mapping of the problem of diffusion in a poten- 
tial onto that of free diffusion with an absorbing region 
for rotational diffusion on a spherical surface. We use an 
attractive Gaussian potential U{6) = — eexp [-{a6f] with 
e = 10 kcal/mol and cr = tt, where the latter corresponds to a 
(half) width of the potential of ^ln(2) /a — 15°, a reasonable 
assumption for a short-range potential. Equating the result- 
ing association rate, Eq. 0, with the rate for free diffusion in 
presence of an absorbing region at 6 = 6o, Eq. we obtain 
Oo ~ 33° for the width of the absorbing region, corresponding 
to an energy drop of AE = U{tt) — U{6o) ~ 0.4 kcal/mol. 



In the case of free diffusion ([/ = 0), Eq. (|SJ becomes 



hiO) - 



(6) 



where we now choose > 0. As before, we equate the 
association rates, Eqs. (0 and ®, to obtain an estimate 
for the width Oq of the absorbing region. 

For a range of parameter values we again find that 
an energy drop of AE = U{n) — U{9n) of the order of 
kT suffices to enter the capture zone (Fig. 0J. We find 
that 9o is insensitive to the choice of e but increases as 
expected with increasing a: the range of the absorbing 
region reflects the range of the potential. 

c. Discussion of the toy model results. The toy mod- 
els have demonstrated that translational and rotational 
diffusional association in a short-range potential (as it 
occurs in protein-protein complex formation when the 
two proteins are close to each other) can be modeled as 
free diffusion in the presence of absorbing regions of suit- 
ably chosen size. The size of the absorbing regions is 
relatively insensitive to the precise shape (that is, the 
functional form) and magnitude of the chosen potential 
function; only the range of the potential must be cho- 
sen roughly right in estimating the angular constraints. 
The energy drop itself required to enter the capture zone 
(binding funnel) is found to be robust towards changes 
in the shape, depth and range of the potential, and can 
therefore be regarded as an essentially universal quantity 
that is largely independent of the particular form of the 
interaction potential used in the mapping problem. 



Since our toy models employ an only one-dimensional 
reaction condition, i.e. a constraint on a single degree 
of freedom, the question arises to what extent the rela- 
tive influence of the potential on the reaction rate would 
change in the case of higher-dimensional reaction condi- 
tions (as used in our subsequent treatment of protein- 
protein interactions where we impose constraints on r, 
9a, B, and 5x)- The results obtained by Zhou (1997) 
suggest that the influence of the interaction potential on 
the association rate constant is more significant for the 
case of two diffusing spheres bearing a circular reactive 
patch on each surface {i.e. where a two-dimensional reac- 
tion condition is used for both spheres) than for the sit- 
uation where one of the spheres is taken to be uniformly 
reactive {i.e. where a two-dimensional reaction condition 
is imposed on one sphere, but an only one-dimensional 
reaction condition is employed for the second sphere). 
Generalizing these findings, we may anticipate that an 
attractive interaction potential will affect reaction rates 
to a larger extent when the number of constrained vari- 
ables in the reaction condition is increased. Since our toy 
models have shown that it suffices to enter the potential 
well by a relatively small amount to be "captured" , we 
can conclude that for the case of a higher-dimensional re- 
action condition as considered in the following, an even 
smaller energy drop will be sufficient to enter the capture 
zone. From the point of view of transition state theory, 
our approach corresponds to identifying the transition 
region and then computing the flux into this region. 



C. Mapping the protein— protein interaction funnel 
from the structure of a protein— protein complex 

Now we shall apply the idea outlined above to an es- 
timate of the angular constraints 9g, 5(j)o, and Sxo, 
needed for the application of our expression for the rate 
constant, Eq. ||2Jl. For this purpose, we have directly 
taken the three-dimensional structures of the considered 
complexes from the Protein Data Bank (PDB). 

First, the side chains of the native complexed struc- 
ture were repacked by minimizing a full-atom energy 
function £ dominated by Lennard- Jones interactions, an 
orientation-dependent hydrogen bond potential, and an 
implicit solvation model (Gray et al. 2003). As with all 
current potential functions for macromolecules, there are 
likely to be considerable inaccuracies in this model, but 
it should be emphasized that the angular constraints and 
rates computed here are relatively insensitive to the de- 
tails of the interactions — the toy examples clearly demon- 
strate that once the binding funnel has been entered 
(which has been found to require only a small drop down 
in energy), the detailed form of the interaction potential 
has only little influence. 

Second, a set of 1,000 alternative structures was gen- 
erated from the native complex by performing random 
small perturbative movements around the native confor- 
mation, and the interaction energy of these structures 
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was evaluated using the same energy function £ as em- 
ployed in the repacking procedure. The energy land- 
scapes defined by these alternative structures exhibit 
clear funnels around the native minimum. 

The toy model calculations show that diffusion in such 
landscapes can be modeled as free diffusion with an effec- 
tive "capture" region several kT into the energy funnels. 
To define the capture energy cutoff £c below which the 
partners are committed to bind, we compute the aver- 
age f av of the energies of the five lowest lying structures 
greater than fO A root mean square deviation (rmsd) 
from the native complex (and hence outside of the na- 
tive energy funnel). Because the energy cutoff cannot be 
determined exactly, we obtain two different estimates of 
the association rate setting to either fav or fav — 5fcT. 
We selected the fO structures with the largest values of 

+ ™ ^"^t of structures with E <Ec and took the 
averages of their values of 9g, S(f>Q, and Sxo to obtain 
estimates for the angular tolerances used in computing 
the association rates. 

An example is shown in Fig. El We see that the funnel- 
like dependence of the energy on the rmsd and on the 
angular deviations is akin to the shape of the attractive 
potentials used in the toy models for translational and 
rotational diffusion discussed above. Furthermore, the 
location of the angular constraints resembles the posi- 
tion of the absorbing regions of the toy models. These 
similarities support our approach of deriving the angular 
constraints 0'^, 0^, 500, and Sxo from the interaction en- 
ergy of perturbed protein complex structures in general, 
and from our method of choosing suitable energy cutoffs 
in particular. 



D. Computation of diffusion-limited association 
rates from structures of protein— protein complexes 

With the angular constraints determined as described 
in the preceding Sec. Ill, we use Eq. ((2Jl in the fully 
diffusion-limited limit (k —>■ oo) to compute protein- 
protein association rates. The effective radii Ra and Rb 
of the spheres representing the proteins are taken to be 
equal to the radius of gyration, Rg — {1/N) J^i di (where 
N is the number of atoms in the protein and di the dis- 
tance of the i-th atom from the geometric center of the 
protein), multiplied by a correction factor of (5/3) to 
obtain the desired result Rg = Rs for the limiting case of 
a homogeneous sphere of radius Rg- The sum Ra + Rb 
is used as the value for the distance R between the cen- 
ters of the two proteins at which reaction is assumed to 
occur, cf. Eqs. Q. The study of our toy model for trans- 
lational diffusion in Sec. II has shown that this serves as 
a good estimate for R since the presence of short-range 
attractive interactions increases the effective reaction ra- 
dius only slightly. The translational and rotational dif- 
fusion constants D = Dj^^"^ + Dg™** and -D™*^, respec- 
tively, are computed from the Stokes-Einstein relations 
DaT = kBT/6TTrjRA,B and D'^'b - kBT/8TTr]R\ g, 
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FIG. 5: Free energy funnels around the native structure. 
The energy £ and rmsd (A), the energy £ and the angular 
deviations O'^ (B) and Sxo (C) are shown for a set of randomly 
perturbed structures of the protein-protein complex IFIN. 
States of lower energy are seen to be associated with smaller 
angles, suggesting that the angles are a reasonable measure 
for the deviation from the correctly complexed structure. The 
two parallel lines represent the two energy cutoffs £c = £a,v 
and £c = £av — 5feT, where fav is the average energy of the five 
lowest energy complexes with a rmsd above 10 A. The vertical 
lines in the plots indicate the resulting angular constraints 6'^ 
and Sxo corresponding to £c — fav (dashed line) and £c = 
£av — 5kT (dotted-dashed line). 
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X 


10* 


2.1 X 10^ 


IGOT 


Gt-a/Gi-Q Chimera 


26.3 


Gt-"'^ 


27.9 


7.5 


4.1 


2.5 


4.1 


5.8 


X 


10* 


1.3 X 10^ 












1.4 


0.6 


0.4 


0.9 


1.2 


X 


10* 


3.4 X 10~* 


IMAH 


AcetylchoUnesterase 


28.3 


Fasciculin-2 


14.2 


7.0 


5.1 


5.2 


0.8 


8.6 


X 


10* 


7.9 












4.0 


3.6 


1.2 


0.2 


3.0 


X 


10* 


7.4 X 10^^ 


ISPB 


SubtiUsin Bpn' 


20.6 


Subtilisin Bpn' 


15.1 


16.9 


17.3 


10.9 


21.1 


2.0 


X 


10"= 


2.6 X 10* 




Prosegment 








12.8 


11.1 


7.2 


13.9 


5.6 


X 


10-"^ 


2.7 X 10^ 


ISTF 


Papain 


20.4 


Papain Inhibitor 


16.4 


8.4 


4.4 


3.9 


1.2 


6.9 


X 


10* 


8.6 








Stefin B 




7.7 


4.0 


3.0 


1.0 


5.4 


X 


10* 


3.8 


ITGS 


Trypsinogen 


20.2 


PSTI 


13.5 


9.1 


4.9 


9.4 


4.2 


1.6 


X 


lO'^ 


1.1 X 10^ 












4.7 


3.5 


4.1 


2.3 


3.6 


X 


10* 


3.5 


2SIC 


SubtiUsin BPN' 


20.9 


Streptomyces Subtilisin 


16.6 


12.1 


18.9 


4.2 


6.6 


5.7 


X 


lO'^ 


1.9 X 10^ 








Inhibitor 




7.0 


13.9 


4.4 


5.5 


2.3 


X 


lO'^ 


3.1 X 10^ 


2TEC 


Thermitase 


21.0 


Eglin-C 


13.9 


8.2 


6.5 


2.5 


9.0 


1.4 


X 


10^ 


8.9 X 10^ 












4.1 


4.1 


2.7 


2.5 


3.2 


X 


10* 


2.6 


3HHR 


Human Growth 


26.1 


Human Growth 


21.2 


13.6 


14.8 


1.9 


3.2 


2.9 


X 


10^ 


3.3 X 10^ 




Hormone 




Hormone Receptor 




10.1 


12.7 


2.0 


0.7 


1.5 


X 


lO'^ 


3.1 X 10^ 


4HTC 


Hirudin 


21.9 


Thrombin 


20.0 


9.3 


7.4 


3.1 


2.8 


9.5 


X 


10* 


5.5 X 10^ 












8.5 


4.9 


3.9 


4.3 


8.2 


X 


10* 


3.9 X 10^ 



TABLE I: Association rates computed from Eq. ((2J in the fully diffusion-controlled limit (k oo) for the set of investigated 
protein-protein complexes. The radii Ra and Rb of each protein in the complex were estimated based on the radius of gyration. 
The angular constraints 6^, 9%, S(f)o, and Sxo were determined as described in the Methods. The geometric rates, shown for 
comparison, are given by kon ~ 4tvDR x d(t>oSxo{i — cos6j^){l 

with 7] = 8.9 X 10""* Ns/m^ (water) and T = 300 K. 

Table I lists the 15 investigated protein-protein inter- 
actions, together with the estimated effective radii Ra 
and Rb, the angular orientational constraints 0^, 0^, 
S(j)o, and 6xoy ^^nd the association rate constants fcon de- 
termined from our theoretical expression, Eq. For 
comparison, we also state the association rates obtained 
from a purely probabilistic model (geometric rates) . 

First of all, it is worth noting that both the angular 
tolerances and the corresponding rate constants are rela- 
tively insensitive (given the approximations involved) to 



- cos e%)/A-R^. 

the particular choice of the energy cutoffs £c = ^av and 
£c = fav — 5fcT. For the protein complexes under study, 
the rates computed from the two energy cutoffs vary in 
average by a factor of 3, and no rates differ by more than 
a factor of 5 for a given complex, thus indicating the 
robustness of our method of estimating these rates. 

We observe that the angular constraints vary signifi- 
cantly among the investigated complexes, which suggests 
that our procedure of estimating these tolerances yields 
indeed characteristic and distinguishable values. 

The association rate constants obtained using these 
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angular constraints range from 10^-10^ s~^ and 

are significantly higher than the corresponding geomet- 
ric rates. While the experimentally determined associa- 
tion rates of many protein-protein complexes are in this 
range, considerably faster rates are also observed, likely 
because of significant long-range interactions neglected 
in our model. 



III. DISCUSSION 

We have presented a simple model for the associa- 
tion of proteins. The molecules are modeled as diffusing 
spheres, no forces are assumed to act between them, and 
the reaction condition is based on an estimate of angu- 
lar constraints on the mutual orientation of the molec- 
ular interfaces based on the assumption of short-range 
guiding forces. This procedure allows for an applica- 
tion of an explicit mathematical expression for the as- 
sociation rate constant that we have previously derived 
(Schlosshauer and Baker 2002). In this paper, we have 
used this method to estimate association rates of a set of 
15 different protein-protein complexes. 

The computed rates all he within 10'*-10'^ M'^ s"\ 
which can thus be taken as the typical diffusion-limited 
protein-protein association rate in the absence of attrac- 
tive interactions, in good agreement with what is exper- 
imentally known for such interactions. This is several 
orders of magnitude higher than the geometric rate that 
had previously been used by various authors. Our result 
therefore shows that typical diffusion-limited association 
rates of proteins where no or or only weak long-range in- 
teractions are present can essentially be explained with a 
model that is solely based on translational and rotational 
diffusion. Experimentally observed significantly higher 
rates typically suggest the presence of electrostatic steer- 
ing forces, whereas much lower rates may indicate a re- 
action that is opposed by free energy barriers and is thus 
not fully diffusion-limited. 

The advantage of our method over the traditional ap- 
proach of BD simulations lies in the fact that our tech- 
nique provides a more physically transparent insight into 
the resulting association rates. The differences in rates 
among protein complexes can directly be traced back to 
the sizes and shapes of the respective reactive zones in 
configurational space, which are determined by mapping 
out the binding funnel in the interaction energy land- 
scape. 

The model completely neglects possible free energy 
barriers due to desolvation and/or side-chain freezing 
during complex formation as well as a possible slowing 
down of diffusion within the binding funnel due to in- 
creased ruggedness of the landscape. Our finding that 
the association rates obtained with the simple diffusional 
model are in the range of those of many protein-protein 
complexes (10^-10^ s~^) suggests that free energy 
barriers and landscape ruggedness do not have a signifi- 
cant impact on the dynamics of protein-protein associa- 



tion. 

Our model provides a zeroth-order estimate of protein- 
protein association rates in the absence of long-range 
interactions. This contrasts with most previous work, 
which has sought to account for changes in association 
rates accompanying sequence changes, rather than the 
absolute association rate. By incorporating long-range 
electrostatic interactions into our diffusional model, it 
should be possible to develop a complete theory of associ- 
ation kinetics that can account for both the sequence de- 
pendence and the absolute magnitude of protein-protein 
association rates. 



Online ressources 

For the interested reader who would like to compute 
protein-protein association rates from our model, we 
have created a web application that allows the user 
to submit the set of reaction conditions parameters, 
Eqs. and then returns the corresponding binding rate 
by evaluating Eq. (|2Jl . The webserver can be accessed via 
http : //tools . bakerlab . org/ ^pprate. 
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